1 Installing and starting the program

Install betAS by typing the following in an R console:

devtools::install_github("marianaferreira/betAS")

After the installation, load betAS and other required packages by typing:

#devtools::load_all() # delete this later, useful to test
library(betAS)
# For visualization
library(ggplot2)
library(plotly) 
library(dplyr)
library(ggpubr)

2 Loading alternative splicing quantification data

betAS accepts files sourced from vast-tools, rMATS, or whippet. These files should contain junction reads information (e.g. results obtained from the vast-tools’ tidy module are not supported).

Based on the value assigned to the tool parameter, the pathTables parameter can take the following input types:

• tool="vast-tools": Should be a path pointing to a *INCLUSION_LEVELS_FULL*.tab file.
• tool="rMATS": Should be a path pointing to a *MATS.JC.txt file.
• tool="whippet": Should be a list of paths (minimum of 2) pointing to *.psi.gz files.

# Example to load data from vast-tools
dataset <- getDataset(pathTables="path/to/dataset/*INCLUSION_LEVELS_FULL*.tab", tool="vast-tools")

# Example to load data from rMATS
dataset <- getDataset(pathTables="path/to/dataset/*MATS.JC.txt", tool="rMATS")

# Example to load data from whippet
dataset <- getDataset(pathTables=list("path/to/sample1/*.psi.gz",
                                      "path/to/sample2/*.psi.gz",
                                      "path/to/sample3/*.psi.gz"), tool="whippet")

For users who may not have files to upload, there’s also an option to work with our built-in dataset. This dataset is a subset of the Deep transcriptional profiling of longitudinal changes during neurogenesis and network maturation in vivo public mouse dataset:

To use the pre-built datasets, simply don’t specify the file and choose the tool you prefer.

# Example to load data from vast-tools
dataset <- getDataset(pathTables=NULL, tool="vast-tools")

# Example to load data from rMATS
dataset <- getDataset(pathTables=NULL, tool="rMATS")

# Example to load data from whippet
dataset <- getDataset(pathTables=NULL, tool="whippet")

3 Filtering datasets

betAS allows to filter the dataset based on alternative splicing event types, PSI value range to consider and minimum number of supporting reads.

For simplicity, in this tutorial we will be using betAS provided dataset from vast-tools to illustrate such features. The getDataset function retrieves the inclusion tables specific to the chosen tool for further analyses, and the getEvents formats inclusion tables into a standard format that is compatible with betAS.

The output from getEvents is an object with 4 attributes:

• PSI: A data frame with information regarding the splicing events and their inclusion levels.
• Qual: A data frame detailing the events and information related to the coverage associated with each event. This representation follows the vast-tools representation of inclusion (inc) and exclusion (exc) reads.
• EventsPerType: A named table summarising the number of events per type.
• Samples: Vector containing the names of the samples considered.

tool <-  "vast-tools"

# Import data
dataset <- getDataset(pathTables=NULL, tool = tool)
# Get data ready for betAS
dataset <- getEvents(dataset, tool = tool)

lapply(dataset,head)

## $PSI
##     GENE        EVENT                     COORD LENGTH                                                      FullCO COMPLEX SRR645824 SRR645826 SRR645828 SRR645830 SRR645872 SRR645875
## 2  Cdc45 MmuEX6061348   chr16:18795095-18795226    132 chr16:18795825+18795831,18795095-18795226,18794860+18794874       S    100.00    100.00    100.00    100.00     100.0    100.00
## 7   Narf MmuEX0030885 chr11:121250293-121250430    138               chr11:121249192,121250293-121250430,121250760       S    100.00    100.00    100.00    100.00     100.0    100.00
## 9  Scmh1 MmuEX0041239  chr4:120524998-120525063     66                chr4:120522649,120524998-120525063,120528156       S    100.00    100.00     92.86    100.00     100.0     94.74
## 11  Tfe3 MmuEX0046712      chrX:7767893-7767997    105                        chrX:7767684,7767893-7767997,7769430       S     93.62     89.36     92.16     87.34      81.4    100.00
## 12 Axin2 MmuEX0007157 chr11:108942935-108943129    195               chr11:108942696,108942935-108943129,108943656       S    100.00    100.00    100.00    100.00     100.0    100.00
## 13 Axin2 MmuEX0007154 chr11:108944510-108944605     96               chr11:108943886,108944510-108944605,108945880       S    100.00    100.00    100.00    100.00     100.0    100.00
##    SRR645877 SRR645879
## 2     100.00    100.00
## 7     100.00    100.00
## 9     100.00     95.35
## 11     94.74     90.24
## 12    100.00    100.00
## 13    100.00    100.00
## 
## $Qual
##     GENE        EVENT                     COORD LENGTH                                                      FullCO COMPLEX                       SRR645824.Q
## 2  Cdc45 MmuEX6061348   chr16:18795095-18795226    132 chr16:18795825+18795831,18795095-18795226,18794860+18794874       S SOK,SOK,92=102=0,OK,S@194.00,0.00
## 7   Narf MmuEX0030885 chr11:121250293-121250430    138               chr11:121249192,121250293-121250430,121250760       S     OK,OK,19=30=0,OK,S@49.00,0.00
## 9  Scmh1 MmuEX0041239  chr4:120524998-120525063     66                chr4:120522649,120524998-120525063,120528156       S  VLOW,VLOW,9=13=0,OK,S@22.00,0.00
## 11  Tfe3 MmuEX0046712      chrX:7767893-7767997    105                        chrX:7767684,7767893-7767997,7769430       S     OK,OK,36=52=3,OK,S@85.19,5.81
## 12 Axin2 MmuEX0007157 chr11:108942935-108943129    195               chr11:108942696,108942935-108943129,108943656       S         N,N,6=8=0,OK,S@14.00,0.00
## 13 Axin2 MmuEX0007154 chr11:108944510-108944605     96               chr11:108943886,108944510-108944605,108945880       S          N,N,5=3=0,OK,S@8.00,0.00
##                           SRR645826.Q                       SRR645828.Q                      SRR645830.Q                      SRR645872.Q                      SRR645875.Q
## 2  SOK,SOK,102=102=0,OK,S@204.00,0.00  SOK,SOK,94=98=0,OK,S@192.00,0.00 SOK,SOK,75=70=0,OK,S@145.00,0.00         N,N,1=2=0,OK,S@3.00,0.00         N,N,2=3=0,OK,S@5.00,0.00
## 7       OK,OK,26=18=0,OK,S@44.00,0.00     OK,OK,26=19=0,OK,S@45.00,0.00    OK,OK,20=20=0,OK,S@40.00,0.00 SOK,SOK,66=68=0,OK,S@134.00,0.00 SOK,SOK,66=54=0,OK,S@120.00,0.00
## 9     LOW,LOW,13=15=0,OK,S@28.00,0.00 VLOW,VLOW,12=14=1,OK,S@25.07,1.93 VLOW,VLOW,12=7=0,OK,S@19.00,0.00  LOW,LOW,12=24=0,OK,S@36.00,0.00  LOW,LOW,17=19=1,OK,S@35.05,1.95
## 11      OK,OK,31=53=5,OK,S@79.53,9.47     OK,OK,44=50=4,OK,S@90.32,7.68    OK,OK,35=34=5,OK,S@64.63,9.37  LOW,LOW,14=21=4,OK,S@31.75,7.25    OK,OK,22=18=0,OK,S@40.00,0.00
## 12           N,N,1=8=0,B2,S@9.00,0.00        N,N,2=13=0,B2,S@15.00,0.00  LOW,LOW,11=16=0,OK,S@27.00,0.00         N,N,4=4=0,OK,S@8.00,0.00         N,N,3=5=0,OK,S@8.00,0.00
## 13           N,N,5=2=0,B1,S@7.00,0.00  VLOW,VLOW,7=10=0,OK,S@17.00,0.00 VLOW,VLOW,12=5=0,B1,S@17.00,0.00        N,N,9=3=0,B1,S@12.00,0.00        N,N,8=4=0,OK,S@12.00,0.00
##                         SRR645877.Q                      SRR645879.Q
## 2          N,N,2=1=0,OK,S@3.00,0.00         N,N,1=2=0,OK,S@3.00,0.00
## 7  SOK,SOK,73=70=0,OK,S@143.00,0.00 SOK,SOK,90=63=0,OK,S@153.00,0.00
## 9     OK,OK,15=24=0,OK,S@39.00,0.00  LOW,LOW,27=14=1,OK,S@40.05,1.95
## 11 VLOW,VLOW,9=27=1,B1,S@35.05,1.95    OK,OK,16=21=2,OK,S@35.19,3.81
## 12         N,N,5=4=0,OK,S@9.00,0.00        N,N,4=8=0,OK,S@12.00,0.00
## 13        N,N,7=5=0,OK,S@12.00,0.00       N,N,10=3=0,B1,S@13.00,0.00
## 
## $EventsPerType
## 
##  Alt3  Alt5   ANN    C1    C2    C3 
## 85993 61115 24750  5699  5179  5896 
## 
## $Samples
## [1] "SRR645824" "SRR645826" "SRR645828" "SRR645830" "SRR645872" "SRR645875"

We will be using betAS metadata for vast-tools built-in example to ease group creation.

data("VT2_metadata_mouse")
(metadata <- as.data.frame(VT2_metadata_mouse))

# to clean-up the environment of unecessary variables, loaded automatically by betAS when using default datasets
rm(VT2_data_mouse)
rm(VT2_metadata_mouse)

dataset_filtered <- filterEvents(dataset, types=c("C1", "C2", "C3", "S", "MIC"), N=10)

cat(paste0("Alternative events: ", nrow(dataset_filtered$PSI)))

## Alternative events: 195636

dataset_filtered <- alternativeEvents(dataset_filtered, minPsi=1, maxPsi=99  )

cat(paste0("Alternative events: ", nrow(dataset_filtered$PSI)))

## Alternative events: 1404

bigPicturePlot <- bigPicturePlot(table = dataset_filtered$PSI)
bigPicturePlot + theme_minimal()

4 Group definition and management from automatically detected or user-defined sample list

To allow for alternative splicing analyses between groups of interest, the user must firstly specify to what groups the samples under analyses belong to.

To check the names of the samples being used, we can access the Samples attribute inside the dataset object.

cat("Sample names: ")

## Sample names:

cat(dataset_filtered$Samples)

## SRR645824 SRR645826 SRR645828 SRR645830 SRR645872 SRR645875 SRR645877 SRR645879

In order to ease group visualization, we can assign a color to each group. In this section we showcase how to assign a group to a given color, based on a manually defined palette. To automatically define color palettes of N colors the user might want to check Custom color palette for group definition section of this tutorial.

samples <- dataset_filtered$Samples 
random_colors <- c("#FF9AA2", "#FFB7B2", "#FFDAC1", "#E2F0CB", "#B5EAD7", "#C7CEEA", "#FBE2FD", "#D9ECFE")

# Automatic group definition
not_grouped   <- samples
checked       <- c()
groups        <- LETTERS[seq(1, length = min(length(names), 26))]
used_groups   <- groups
groupList <- list()

while((length(checked) < length(names)) & length(used_groups) <= 26){
  
  groupNames  <- agrep(pattern = not_grouped[1], x = not_grouped, value = TRUE)
  checked     <- c(checked, groupNames)
  not_grouped <- not_grouped[-c(match(groupNames, not_grouped))]
  
  # Assign new group
  currentNames <- names(groupList)
  groupList[[length(groupList)+1]] <- list(name = used_groups[1],
                                                   samples = groupNames,
                                                   color = random_colors[1])
  names(groupList) <- make.unique(c(currentNames, used_groups[1]))
  
  random_colors <- random_colors[-1]
  used_groups   <- used_groups[-1]
  
}

groupList

## $A
## $A$name
## [1] "A"
## 
## $A$samples
## [1] "SRR645824" "SRR645826" "SRR645828"
## 
## $A$color
## [1] "#FF9AA2"
## 
## 
## $B
## $B$name
## [1] "B"
## 
## $B$samples
## [1] "SRR645830"
## 
## $B$color
## [1] "#FFB7B2"
## 
## 
## $C
## $C$name
## [1] "C"
## 
## $C$samples
## [1] "SRR645872" "SRR645875" "SRR645877" "SRR645879"
## 
## $C$color
## [1] "#FFDAC1"

# Visualize colors being used
slices <- rep(1, length(groupList))  # Equal-sized slices for each color
# Display the pie chart with colors
pie(slices, col = random_colors[1:length(groupList)], border = "black", labels=names(groupList), main = "Color palette for group definition")

# Define variable of the metadata table to be used as a grouping variable
groupingVariable <- "CellType"
# to get the unique values of the grouping variable. In this example, will be Neuron and ESC
groups <- unique(metadata[,groupingVariable])
# Define vector of sample names based on the example metadata 
samples <- metadata$Run
# Define colors for the two groups
random_colors <- c("#FF9AA2", "#FFB7B2")

groupList <- list()

for(i in 1:length(groups)){

  groupNames <- samples[which(metadata[,groupingVariable] == groups[i])]

  # Assign new group
  currentNames <- names(groupList)
  groupList[[length(groupList)+1]] <- list(name = groups[i],
                                           samples = groupNames,
                                           color = random_colors[i])
  names(groupList) <- make.unique(c(currentNames, groups[i]))

}

groupList

## $ESC
## $ESC$name
## [1] "ESC"
## 
## $ESC$samples
## [1] "SRR645824" "SRR645826" "SRR645828" "SRR645830"
## 
## $ESC$color
## [1] "#FF9AA2"
## 
## 
## $Neuron
## $Neuron$name
## [1] "Neuron"
## 
## $Neuron$samples
## [1] "SRR645872" "SRR645875" "SRR645877" "SRR645879"
## 
## $Neuron$color
## [1] "#FFB7B2"

# Visualize colors being used
slices <- rep(1, length(groups))  # Equal-sized slices for each color
# Display the pie chart with colors
pie(slices, col = random_colors[1:length(groups)], border = "black", labels=groups, main = "Color palette for group definition")

5 Differential alternative splicing between two groups

betAS approach for alternative splicing analyses takes into account two metrics for each alternative splicing event under analyses: the magnitude of the effect (𝚫PSI = PSI_betAS(group A) - PSI_betAS(group B)) and the statistical significance of such changes. For the latter, betAS introduces three approaches: probability of differential splicing (Pdiff), F-statistic or false discovery rate (FDR).

• Pdiff: This approach takes the two sets of random points per condition and calculates, for each AS event’s estimated ∆PSI, the proportion of differences between these that are greater than zero, which has the same interpretation as asking what proportion of beta distribution-emitted values for one condition are higher than those emitted for the other, thus reflecting the probability of differential AS of PSI_betAS(group A) being greater than PSI_betAS(group B).

• F-statistic: betAS also enables an ANOVA-like analysis of variance, comparing inter- and intra-group variabilities. For each event, within is considered the set of differences between each pair of samples that are part of the same group and between the set of differences between each pair of groups. The ratio of the median absolute values of between and within therefore provides an “F-like” statistic. This metric provides a compromise between the effect size of AS differences and their significance.

• FDR: Random generation of points from a beta distribution is used to estimate the null distribution’s PSI and its precision. Then, one point is randomly selected from each of the sample’s null distribution and, keeping the samples’ group assignment (i.e., which samples belong to each group), the ∆PSI between groups under the null hypothesis is calculated. The process is repeated many times (10 000 by default) and the FDR is the proportion of ∆PSI random simulations that are larger than (i.e., more extreme) or equal to the empirical ∆PSI.

# Define groups
groupA    <- "Neuron"
groupB    <- "ESC"

# Define samples inside each group
samplesA    <- groupList[[groupA]]$samples
samplesB    <- groupList[[groupB]]$samples

# Convert samples into indexes
colsGroupA    <- convertCols(dataset_filtered$PSI, samplesA)
colsGroupB    <- convertCols(dataset_filtered$PSI, samplesB)


cat("samplesA:")

## samplesA:

cat(samplesA) 

## SRR645872 SRR645875 SRR645877 SRR645879

cat("\nsamplesB:")

## 
## samplesB:

cat(samplesB)

## SRR645824 SRR645826 SRR645828 SRR645830

cat("\n")

cat("colsGroupA: \n")

## colsGroupA:

colsGroupA

## SRR645872 SRR645875 SRR645877 SRR645879 
##        11        12        13        14

cat("colsGroupB:  \n")

## colsGroupB:

colsGroupB

## SRR645824 SRR645826 SRR645828 SRR645830 
##         7         8         9        10

volcanoTable_Pdiff <- prepareTableVolcano(psitable = dataset_filtered$PSI,
                                    qualtable = dataset_filtered$Qual,
                                    npoints = 500,
                                    colsA = colsGroupA,
                                    colsB = colsGroupB,
                                    labA = groupA,
                                    labB = groupB,
                                    basalColor = "#89C0AE",
                                    interestColor = "#E69A9C",
                                    maxDevTable = maxDevSimulationN100)

volcano_Pdiff <- plotVolcano(betasTable = volcanoTable_Pdiff,
                            labA = groupA,
                            labB = groupB,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") 

volcano_Pdiff

## Warning: ggrepel: 743 unlabeled data points (too many overlaps). Consider increasing max.overlaps

# change size points manually
volcano_Pdiff$layers[[1]]$aes_params$size <- 1
volcano_Pdiff$layers[[2]]$aes_params$size <- 1

# Plot interactively
plotly_volcano_Pdiff <- ggplotly(volcano_Pdiff, width = 700, height = 500) %>%
                          layout(
                            font = list(size = 10),
                            xaxis = list(
                              title = list(font = list(size = 14)),  # Adjust as necessary
                              tickfont = list(size = 10)  # Adjust the tick font size here
                            ),
                            yaxis = list(
                              title = list(font = list(size = 14)),  # Adjust as necessary
                              tickfont = list(size = 10)  # Adjust the tick font size here
                            )
                          )

## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

hover_text  <- paste(
  "Event: ", volcanoTable_Pdiff$EVENT,
  "<br>Pdiff: ", round(volcanoTable_Pdiff$Pdiff,3),
  "<br>Deltapsi: ", round(volcanoTable_Pdiff$deltapsi,3),
  sep = ""
)

 

plotly_volcano_Pdiff$x$data[[1]]$text <- hover_text # all points
plotly_volcano_Pdiff$x$data[[2]]$text <- NULL 
  
plotly_volcano_Pdiff

head(volcanoTable_Pdiff[,c("GENE","EVENT","COORD","Pdiff","deltapsi")])

volcanoTable_Fstat <- prepareTableVolcanoFstat(psitable = dataset_filtered$PSI,
                                    qualtable = dataset_filtered$Qual,
                                    npoints = 500,
                                    colsA = colsGroupA,
                                    colsB = colsGroupB,
                                    labA = groupA,
                                    labB = groupB,
                                    basalColor = "#89C0AE",
                                    interestColor = "#E69A9C",
                                    maxDevTable = maxDevSimulationN100)

volcano_Fstat <- plotVolcanoFstat(betasTable = volcanoTable_Fstat,
                            labA = groupA,
                            labB = groupB,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") 

volcano_Fstat

## Warning: ggrepel: 719 unlabeled data points (too many overlaps). Consider increasing max.overlaps

head(volcanoTable_Fstat[,c("GENE","EVENT","COORD","Fstat","deltapsi")])

volcanoTable_FDR <- prepareTableVolcanoFDR(psitable = dataset_filtered$PSI,
                                    qualtable = dataset_filtered$Qual,
                                    npoints = 500,
                                    colsA = colsGroupA,
                                    colsB = colsGroupB,
                                    labA = groupA,
                                    labB = groupB,
                                    basalColor = "#89C0AE",
                                    interestColor = "#E69A9C",
                                    maxDevTable = maxDevSimulationN100,
                                    nsim = 100) 
        
volcano_FDR <- plotVolcanoFDR(betasTable = volcanoTable_FDR,
                            labA = groupA,
                            labB = groupB,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") 

volcano_FDR

## Warning: ggrepel: 734 unlabeled data points (too many overlaps). Consider increasing max.overlaps

head(volcanoTable_FDR[,c("GENE","EVENT","COORD","FDR","deltapsi")])

eventID <- "MmuEX0003638"

# Auxiliary data to add increment based on the event coverage to avoid dividing by zero 
data("maxDevSimulationN100")

tdensities <- plotIndividualDensitiesList(eventID = eventID,
                                          npoints = 500,
                                          psitable = dataset$PSI,
                                          qualtable = dataset$Qual,
                                          groupList = groupList,
                                          maxDevTable = maxDevSimulationN100)

tdensities + theme_minimal() + ggtitle(eventID)

## Picking joint bandwidth of 0.0033

## Picking joint bandwidth of 0.0135

subset(volcanoTable_Pdiff[,c("GENE","EVENT","COORD","Pdiff","deltapsi")], EVENT==eventID)

subset(volcanoTable_Fstat[,c("GENE","EVENT","COORD","Fstat","deltapsi")], EVENT==eventID)

subset(volcanoTable_FDR[,c("GENE","EVENT","COORD","FDR","deltapsi")], EVENT==eventID)

plotPdiff <- prepareTableEvent(eventID = eventID,
                               psitable = dataset$PSI,
                               qualtable = dataset$Qual,
                               npoints = 500,
                               colsA = colsGroupA,
                               colsB = colsGroupB,
                               labA = groupA,
                               labB = groupB,
                               basalColor = "#89C0AE",
                               interestColor = "#E69A9C",
                               maxDevTable = maxDevSimulationN100,
                               nsim = 1000) %>% 
  plotPDiffFromEventObjList()

plotFstat <- prepareTableEvent(eventID = eventID,
                               psitable = dataset$PSI,
                               qualtable = dataset$Qual,
                               npoints = 500,
                               colsA = colsGroupA,
                               colsB = colsGroupB,
                               labA = groupA,
                               labB = groupB,
                               basalColor = "#89C0AE",
                               interestColor = "#E69A9C",
                               maxDevTable = maxDevSimulationN100,
                               nsim = 1000) %>% 
  plotFstatFromEventObjList()

plotFDR <- prepareTableEvent(eventID = eventID,
                               psitable = dataset$PSI,
                               qualtable = dataset$Qual,
                               npoints = 500,
                               colsA = colsGroupA,
                               colsB = colsGroupB,
                               labA = groupA,
                               labB = groupB,
                               basalColor = "#89C0AE",
                               interestColor = "#E69A9C",
                               maxDevTable = maxDevSimulationN100,
                               nsim = 1000) %>% 
  plotFDRFromEventObjList()

ggarrange(plotPdiff,plotFstat,plotFDR, ncol=3)

6 Differential alternative splicing between multiple groups

The differential AS approach implemented by betAS can be applied to multiple (i.e., more than two) groups in a novel ANOVA-inspired way that extends the Pdiff definition to the comparison of the differences between samples belonging to different biological conditions to those found between replicates. To illustrate this, we applied multiple-group betAS to the analysis of AS differences in a subset of human transcriptomes of forebrain, hindbrain, heart, kidney, liver and testis. This dataset is a subset of the Human RNA-seq time-series of the development of seven major organs public human dataset.

# load data and prepare it for betAS
data("VT1_data_human")
dataset_multgroups <- getEvents(VT1_data_human, tool = tool)

# load metadata
data("VT1_metadata_human")
metadata_multgroups <- as.data.frame(VT1_metadata_human)
metadata_multgroups

# Define variable of the metadata table to be used as a grouping variable
groupingVariable <- "organism_part"
# to get the unique values of the grouping variable. In this example, will be Neuron and ESC
groups <- unique(metadata_multgroups[,groupingVariable])
# Define vector of sample names based on the example metadata 
samples <- metadata_multgroups$Run
# Define colors for the two groups
random_colors <- c("#FF9AA2", "#FFB7B2", "#FFDAC1", "#E2F0CB", "#B5EAD7")

groupList <- list()

for(i in 1:length(groups)){

  groupNames <- samples[which(metadata_multgroups[,groupingVariable] == groups[i])]

  # Assign new group
  currentNames <- names(groupList)
  groupList[[length(groupList)+1]] <- list(name = groups[i],
                                           samples = groupNames,
                                           color = random_colors[i])
  names(groupList) <- make.unique(c(currentNames, groups[i]))

}

groupList

## $forebrain
## $forebrain$name
## [1] "forebrain"
## 
## $forebrain$samples
##  [1] "ERR2598266" "ERR2598267" "ERR2598268" "ERR2598269" "ERR2598280" "ERR2598281" "ERR2598305" "ERR2598306" "ERR2598346" "ERR2598347" "ERR2598348" "ERR2598349" "ERR2598350"
## 
## $forebrain$color
## [1] "#FF9AA2"
## 
## 
## $heart
## $heart$name
## [1] "heart"
## 
## $heart$samples
## [1] "ERR2598270" "ERR2598307" "ERR2598351"
## 
## $heart$color
## [1] "#FFB7B2"
## 
## 
## $hindbrain
## $hindbrain$name
## [1] "hindbrain"
## 
## $hindbrain$samples
##  [1] "ERR2598271" "ERR2598272" "ERR2598273" "ERR2598274" "ERR2598282" "ERR2598283" "ERR2598308" "ERR2598309" "ERR2598352" "ERR2598353" "ERR2598354" "ERR2598355" "ERR2598356"
## 
## $hindbrain$color
## [1] "#FFDAC1"
## 
## 
## $liver
## $liver$name
## [1] "liver"
## 
## $liver$samples
## [1] "ERR2598275" "ERR2598284" "ERR2598285" "ERR2598310" "ERR2598357" "ERR2598358" "ERR2598359"
## 
## $liver$color
## [1] "#E2F0CB"
## 
## 
## $testis
## $testis$name
## [1] "testis"
## 
## $testis$samples
##  [1] "ERR2598276" "ERR2598277" "ERR2598278" "ERR2598279" "ERR2598286" "ERR2598311" "ERR2598360" "ERR2598361" "ERR2598362" "ERR2598363"
## 
## $testis$color
## [1] "#B5EAD7"

# Visualize colors being used
slices <- rep(1, length(groups))  # Equal-sized slices for each color
# Display the pie chart with colors
pie(slices, col = random_colors[1:length(groups)], border = "black", labels=groups, main = "Color palette for group definition")

7 Differential alternative splicing across multiple groups/conditions

8 Supplementary Information

#install.packages("grDevices")
library(grDevices)

generate_pastel_colors <- function(N) {
  pastel_colors <- vector("list", length = N)
  
  for (i in 1:N) {
    # Generate random values for hue, chroma, and luminance
    hue <- runif(1, min = 0, max = 360)  # Adjust the range for a different initial hue
    chroma <- runif(1, min = 50, max = 80)  # Adjust these values for your desired pastel range
    luminance <- runif(1, min = 80, max = 90)  # Adjust these values for your desired pastel range
    
    # Create a pastel color using hcl() with the random values
    pastel_color <- hcl(hue, chroma, luminance)
    
    pastel_colors[[i]] <- pastel_color
  }
  
  return(pastel_colors)
}


# Generate 10 random pastel colors
N <- 10
random_pastel_colors <- generate_pastel_colors(N)

cat("HEX codes of generated pastel colors: ")

## HEX codes of generated pastel colors:

cat(unlist(random_pastel_colors))

## #FFC2DF #B5F07A #B1EB9D #BAE696 #85E6A7 #FFCBA1 #FFB4DA #8FD2FF #55EABF #8FCFFF

# Prepare data for the pie chart
slices <- rep(1, N)  # Equal-sized slices for each color
# Display the pie chart with pastel colors
pie(slices, col = unlist(random_pastel_colors), border = "black", labels=unlist(random_pastel_colors), main = "Automatically generated color palette for group definition")